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Abstract. We simulate the propagation of a planar crack in a quasi-two dimensional fuse model, confining 
the crack between two horizontal plates. We investigate the effect on the roughness of microcrack nucleation 
ahead of the main crack and study the structure of the damage zone. The two dimensional geometry 
introduces a characteristic length in the problem, limiting the crack roughness. The damage ahead of the 
crack does not appear to change the scaling properties of the model, which are well described by gradient 
percolation. 

PACS. 6 2.20.Fe, 62.20.Mk, 64.60.Lx 



1 Introduction 



Understanding the mechanisms of crack propagation is 
an important issue in mechanics, with potential appli- 
cation to geophysics and material science. Experiments 
have shown that in several materials under different load- 
ing conditions, the crack front tends to roughen and can 
often be described by self-affine scaling In particular, 
the out of plane roughness exponent is found displays uni- 
versal values for a wide variety of materials || . Interesting 
experiments have been recently performed on PMMA and 
the in plane roughness of a planar crack was observed to 
scale with an exponent ( — 0.63 ± 0.03 |J. 

While the experimental characterization of crack rough- 
ness is quite advanced and the numerical results very ac- 
curate, theoretical understanding and numerical models 
are still unsatisfactory. The simplest theoretical approach 
to the problem identifies the crack front with a deformable 
line pushed by the external stress through a random tough- 
ness landscape. The deviations of the crack front from a 
flat line are opposed by the elastic stress field, through the 
stress intensity factor 0]. In certain conditions, the prob- 
lem can be directly related to models and theories of in- 
terface depinning in random media and the roughness ex- 
ponent computed by numerical simulations and renormal- 
ization group calculations |^||. Unfortunately, the agree- 
ment within this theoretical approach and experiments are 
quite poor. For the out of plane roughness the theory pre- 
dicts only a logarithmic roughness in mode I |Q , while the 
experimental results give (± — 0.5 at small length scale 
and C± = 0.8 at larger length scales ||. For planar cracks, 
simulations predict £ = 0.35 || and RG gives £ = 1/3 ||, 
both quite far from the experimental result. The inclu- 



sion of more details in the model, such as elastodynamic 
effects, does not lead to better results ||. 

A different approach to crack propagation in disor- 
dered media, considers the problem from the point of view 
of lattice models . The elastic medium is replaced by 
a network of bonds obeying discretized equations until 
the stress reached a failure threshold. The disorder in the 
medium can be simulated by a distribution of thresholds 
or by bond dilution. Models of this kind have been widely 
used in the past to investigate several features of fracture 
of disordered media, such as the failure stress 11 12 , 13 , 14 1 
fractal properties p0|Jlr| and avalanches U[ 16,17,18,191 



The out of plane roughness exponent has been simulated 
in two dimension, resulting in £j_ ~ 0. 7 pO|,pT| , and three 
dimensions where (± ~ 0.4 - 0.5 ^[^jTThe last re- 
sult is in good agreement with experimental results, if we 
identify the small length scales with the quasistatic regime 
used in simulations. The advantage of lattice models over 
interface models is that the former allow for nucleation 
of microcracks ahead of the main crack. While it is well 
known experimentally that microcracks do nucleate, their 
effect on the roughness exponent has never been studied. 

In this paper we present numerical simulations of a pla- 
nar crack using the random fuse model JlT[ . We employ 
a quasi two-dimensional geometry, considering two hor- 
izontal plates separated by a network of vertical bonds. 
A similar setup was used in a spring model p5[ , but the 
roughness was studied only in the high velocity regime 
for crack motion. The experiments of Ref. || where in- 
stead performed at low velocity so that a quasistatic model 
seems more appropriate. 

We find that the two dimensional geometry introduces 
a characteristic length limiting the crack roughness. In 
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addition, crack nucleation does not appear to change in 
a qualitative way the behavior of the system. For length 
scales smaller than the characteristic length, the crack is 
not self affine, but possibly self-similar. We study the dam- 
age zone close to the crack and find that several of its 
features can be described by gradient percolation p7|. 



2 Model 

In the random fuse model, each bond of the lattice repre- 
sents a fuse, that will burn when its current overcome a 
threshold ]ll| , p^ , p^|JT^ ] . The currents flowing in the lattice 
are obtained solving the Kirchhoff equation with appro- 
priate boundary conditions p8| . In this paper, we consider 
two horizontal tilted square lattices of resistors connected 
by vertical fuses (see Fig. [I]). The conductivity of the hor- 
izontal resistors is chosen to be unity, while the vertical 
fuses have conductivity a. A voltage drop AV is imposed 
between the first horizontal rows of the plates. To simu- 
late the propagation of a planar crack, we allow for fail- 
ures of vertical bonds only and assign to each of them a 
random threshold jf, uniformly distributed in the interval 
[1 : 2). When the current in a bond i overcomes the ran- 
dom threshold, the bond is removed from the lattice and 
the currents in the lattice are recomputed, until all the 
currents are below the threshold. The voltage drop is thus 
increased until the weakest bond reaches the threshold. 

The quasistatic dynamics we are using should cor- 
respond to the small constant displacement rate at the 
boundary of the crack used in experiments |3) . In order to 
avoid spurious boundary effect, we start with a preexist- 
ing crack occupying the first half of the lattice (see Fig. 
and employ periodic boundary conditions in the direction 
parallel to the crack. In addition, once an entire row of 
fuses has failed, we shift the lattice backwards one step in 
the direction perpendicular to the crack, to keep the crack 
always in the middle of the lattice. 

Before discussing the numerical results, we present some 
analytical considerations which will guide the simulations. 



3 Characteristic length 

Here we investigate the model introduced in the preceding 
section in some particular configuration. We first analyze 
the case of a perfectly straight planar crack and study the 
current decay in front of it. In this condition, the system 
is symmetric in the direction parallel to the crack and we 
can thus reduce it to one dimension. 

We consider an infinite ladder composed of vertical 
bonds of resistance r = 1/a connected by unitary hori- 
zontal resistances. Since the ladder is infinite, we can add 
one additional step without changing the end to end re- 
sistance R: 



R = 2 + l/(l/r+l/R) = 2 + 



rR 
r + R' 



Solving Eq. (|l]) we obtain the total resistance 
R = VT+2r + l. 



(1) 



(2) 



The fraction of current j flowing through the first ladder 
step is such that rj = R(l — j), which implies 



R 



VT+2r+ 1 



vT+27- 1 



R 



r + l + Vl + 27 



(3) 



The current flowing in the second ladder step is then (1 — 
j)j and similarly in the nth step it is given by j n = (1 — 
j)" j, thus scaling as j n cx exp(— n/£) where 



-i 



-l 



log(l - j) log(l - tr(l - y/l-2/a)) 



(4) 



Thus the current in front of the crack decays exponen- 
tially with a characteristic length £ ~ 1/ \2a, for a <C 1. A 
similar result could have been anticipated from the struc- 
ture of the Kirchhoff equations reading as 



5>* 



i-\-nn 



Vi) + a{V! - Vi) = o 



(5) 



where the sum runs over the nearest neighbors of node 
i and V[ is the voltage of the corresponding node in the 
opposite plate. Due to symmetry we can chose V( = —Vi 
and solve the equations only for one of the plates. Eq. (||) 
represents a discretization of a Laplace equation with a 
"mass term" 



v'v-c v = °> 



(6) 



where £ = 1/ \/2a. 

The continuum limit can be used to understand how 
current is transfered after a single failure. We define G{x— 
x') as the difference in the currents in x before and after 
a bond in x' has failed. The function G is analogous of 
the "stress Green function" used in interface models for 
cracks propagation In fact, the equation of motion 

in these models is written as 



dh(x, t) 
dt 



F+ I dyG(x-y)(h(y,t)-h(x,t))+r](x,h), 

(7) 

where h indicate the position of the crack, F is propor- 
tional to the external stress, rj to the random toughness 
of the material and for a planar crack in three dimensions 
G(x) ~ l/|a;| 2 . A renormalization group analysis shows 
that the roughness of the interface crucially depends on 

\ for 



the decay of G 



If G decays slower than 



x — ► oo the interface is not rough on large length scales. 

In order to compute the function G, we solve Eq. (^) 
with the appropriate boundary conditions. Note that by 
definition G(x) is proportional to the differences in the 
voltages in x before and after removing a fuse. Since Eq. (|J) 
is linear, the difference of the voltages still satisfies the 
equation in all points except x = 0. This condition can 
also be expressed in terms of the current J = dV/ dy p6[ , 
which should be continuous everywhere apart from x = 0. 

Let's consider a planar crack along the x direction and 
identify two domains: 

1) the domain where fuses are present (y > 0) labeled A 

2) the domain where all fuses are burnt out (y < 0) labeled 
B. Thus the equation to solve in domain A is 



w 2 v = c v 



(8) 
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and in B V 2 V = 0. 

Taking the Fourier transform along x, and calling k 
the conjugate variable to x, we can write in domain A 



d 2 v = (k 2 + c 2 )v. 



(9) 



Integrating the equation, setting V — > at infinity, we 
obtain 

V(Ky) = V(k,0)exp(-y/£) (10) 



where l/£ — \Jk 2 + £ 2 . A similar calculation allows to 
obtain V(k, y) in domain B. 

The currents normal to the crack in the two domains 
are given by 



and 



J A {k) = -Vk 2 +^- 2 V(k,0) 



J B {k) = -\k\V(k,0), 



(11) 



(12) 



where V is the same for the two domains. If one bond is 
removed at x — along the interface, the continuity of 
the current implies J a + Jb ~ 5(x) and in Fourier space 



Mk) + J B (k) ~ 1 



and hence 



V 



-1 



1*1 + Vf^+f 



(13) 



(14) 



The Fourier transform of the function G is simply pro- 
portional to V and therefore at short distances k£ 1, 
V oc l/\k\, or G oc log(a;), while at long distances fc£ <C 1, 
G oc exp(— r/£). 

We test the asymptotic behavior predicted above by 
estimating G from numerical simulations. The results for 
a lattice of size L = 128 are in good agreement with the 
analytical predictions as shown in Fig. [3| It is interesting 
to remark that the roughness of the crack is limited by 
£, but even in the limit £ — > oo we do not expect a self 
affine crack, since G decays slower than l/\x\. In the next 
section, we will show numerically that damage nucleation 
does not alter this conclusion. 



4 Crack roughness: simulations 

In order to analyze the effect of crack nucleation ahead 
of the main crack, we first simulate the model confin- 
ing the ruptures on the crack surface. In this way, our 
model reduces to a connected interface moving in a ran- 
dom medium with an effective stiffness given by the solu- 
tion of the Kirchhoff equations. The results are then com- 
pared with simulations of the unrestricted model, where 
ruptures can occur everywhere in the lattice. In both cases 
the crack width increases with time up to a crossover time 
at which it saturates. In Fig. |^ we compare the damage 
structure in the saturated regime for the two growth rules. 
The height- height correlation function C{x) = ((h(x) — 
h(0)) 2 ), where the average over different realizations of 
the disorder, is shown in Fig. o. From these figures it is 



apparent that the structure of the crack is similar in the 
two cases. The only difference lies in the higher saturation 
width that is observed when microcracks are allowed to 
nucleate ahead of the main crack. 

Next, we analyze the behavior of the crack as a func- 
tion of a which should set the value of the characteristic 
length to £ ~ \j\f2o. In this study we restrict our at- 
tention to the general model with crack nucleation. We 
compute the global width W = ({h 2 ) — (/i) 2 ) 1 / 2 , averag- 
ing over several realization of the disorder (typically 10), 
as a function of time for different values of £. Fig. ^ shows 
that W increases linearly in time until saturation. The 
global width in the saturated regime scales as with 
£ ~ 0.75, as shown in Fig. ^. Due to the limited scaling 
range, we could not obtain a more reliable estimate of the 
exponent value. 



5 Mean-field approach 

The long-range nature of the Green function suggests that 
a mean-field approach could be suitable. We outline here 
the spirit of such an approach, through the determination 
of the density of burnt fuses ahead of the crack front. First, 
we note that the mean profile is expected to be transla- 
tional invariant along the x axis and thus the problem 
reduces to a one dimensional geometry. As argued ear- 
lier, the tension V(y) should obey the following differential 
equation, in the continuum limit: 



d 2 V 
dy 2 



2a(y)V(y) 



(15) 



where <j{y) is the (rc-averaged) conductivity at position y. 
The latter can be written as (1 — D(y))ao where D{y) is 
the "damage", i.e. fraction of burnt fuses. This fraction is 
a known function of the current density in the mean-field 
approach. Namely, the vertical current going through in- 
tact fuses is j(y) = 2a V(y) — V{y)/^ 2 . It is remarkable 
that D(y) drops out of this equation: the damage reduces 
the current density flowing at a given position by a fac- 
tor (1 — D(y)) as compared to the intact state, but the 
same current density flows through a reduced number of 
intact fuses, and is thus multiplied by 1/(1 — D(y)). In 
conclusions the two factors cancel out. 

The proportion of fuses which may support the current 
without burning is given by the cumulative distribution 
of threshold currents: P(j) = Jj p(j')dj' , which implies 
D{y) = 1 — P(j(y)). The voltage profile along the y axis 
is thus given by 



.d 2 v 

dy 2 



P(V(y)/e)V(y) 



(16) 



This equation can be rewritten in terms of the rescaled 
coordinate s — y/£ and current j — V/£ 2 and in our case, 
since P(j) — 2 — j for 1 < j < 2, we obtain 

j"(s)=j(s)(2-j(s)). (17) 
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Notice that Eq. [l?] is valid only for 1 < j < 2, while 
for j < 1 the equation becomes j" — j and for j > 2 
we have j" — 0. At infinity, j < 1 thus the current is 
given by j = e~( s ~ s °\ for s > s , so that the bound- 
ary condition for Eq. 17 at s = sq is j(so) = —j'(so) = 
1 With these boundary conditions, Eq. ( J17| ) can not be 
solved explicitely so we resort to numerical integration. 
From the solution of Eq. ( |l7|) we obtain the damage pro- 
file and compare it to numerical simulations (see Fig. ||) 
The remarkable agreement between the mean-field solu- 
tion and simulations, with a damage profile which is a sin- 
gle function of y/£, implies that the fracture front should 
be given by gradient percolation (in fact the gradient is 
non- linear) [£7] . From this observation we can extract the 
scaling of the front with the gradient g (here g oc l/£) as 



W cx g 



£v/(i+v) w ;h ere v j s the percolation cor- 



relation length critical exponent v = 4/3, or W oc £ ' 57 , 
reasonably consistent with our data. 



6 Conclusions 

In this paper, we have studied the propagation of planar 
cracks in the random fuse model. This model allows to 
investigate the effect on the crack front roughness of the 
microcracks nucleating ahead of the main crack. The study 
was restricted to a quasi two dimensional geometry and 
could apply to cases in which the material is very thin in 
the direction perpendicular to the crack plane [E9[ . 

In two dimensions, the geometry of the lattice induces 
a characteristic length £ limiting the roughness and mi- 
crocrack nucleation does not appear to be relevant. In ad- 
dition, for length scales smaller than £ the Green function 
decays very slowly, suggesting the validity of a mean-field 
approach. We study the problem numerically, computing 
the scaling of the crack width with time and £, and ana- 
lyze the damage ahead of the crack. The results suggest 
an interpretation in terms of gradient percolation [|27| , as 
it is also indicated by mean-field theory. The limited range 
of system sizes accessible to simulations does not allow for 
a definite confirmation of these results. 

The present analysis does not resolve the issue of the 
origin of the value of the roughness exponent for planar 
cracks in heterogeneous media. While microcrack nucle- 
ation is irrelevant in the present context, three dimen- 
sional simulations are needed to understand whether this 
is true in general. In principle, one could still expect that 
microcrack nucleation in three dimensions would change 
the exponent of the interface model (£ = 1/3), but the 
present results do not lead to such a conclusion. 
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Fig. 1. The geometry of the model. The horizontal bonds have 
unitary conductivity while the vertical bond have conductivity 
a. A planar crack is present at the center of the system and a 
voltage drop is applied at the boundaries. 
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Fig. 2. The one dimensional ladder model used to compute the 
characteristic length. The vertical bonds have resistance r the 
horizontal have unitary resistance. The end to end resistance 
of the ladder is R. 
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Fig. 3. The current transfer function in log-log plot (top), to 
show the logarithmic behavior at short length scales. In the 
linear-log plot (bottom) the function G has been rescaled with 
£ and the collapse is in agreement with the analytical solution. 
The deviation from a pure exponential exp(—x/£), plotted as 
a dashed line, is due to the periodic boundary conditions em- 
ployed in the simulation. 



6 



Stefano Zapperi et al.: Planar cracks in the fuse model 



100 



80 



(a) 



ft - _ sr-. .* i— s" 
rv» - - —v s 'ja -a 



?S 60 



40 



20 



°JS?° JJig rest fart: ;«« 



50 



100 



4.9 



4.7 - 



4.5 



4.3 - 



4.1 



3.9 



10 

X 



100 



100 



?s 60 




Fig. 4. The crack profiles in the steady state, when only a 
connected crack is allowed to grow (a) and when damage is 
not restricted to occur only close to the crack (b). The crack 
interface is shown in dark. 



Fig. 5. Comparison between the height correlation function of 
a connected crack (bottom curve) and that of a unrestricted 
crack (upper curve). 
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Fig. 6. The global width as a function of time for different 
values of a and £ = 1/ v2o"- 
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Fig. 7. The saturated global width as a function of £. 
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Fig. 8. The average concentration D(y) of burnt fuses as a 
function of the reduced distance y/£ from the crack for different 
values of a. The numerical results are in excellent agreement 
with the mean-field solution. 



